Correlation potentials for molecular bond dissociation within the self-consistent 

random phase approximation 
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Self-consistent correlation potentials for H2 and LiH for various inter-atomic separations are ob- 
tained within the random phase approximation (RPA) of density functional theory. The RPA 
correlation potential shows a peak at the bond midpoint, which is an exact feature of the true cor- 
relation potential, but lacks another exact feature: the step important to preserve integer charge on 
the atomic fragments in the dissociation limit. An analysis of the RPA energy functional in terms 
of fractional charge is given which confirms these observations. We find that the RPA misses the 
derivative discontinuity at odd integer particle numbers but explicitly eliminates the fractional spin 
error in the exact-exchange functional. The latter finding explains the accurate total energy in the 
dissociation limit. 



I. INTRODUCTION 

The random phase approximation (RPA) in Kohn- 
Sham (KS) density functional theory (DFT) has in re- 
cent years received considerable attention in quantum 
chemistry^"'^ and material science. ''"^^ The RPA incorpo- 
rates exchange effects exactly and the correlation energy 
is treated non-perturbatively by summing a subset of po- 
larization diagrams to infinite order. ^^'^'^ Furthermore, 
the RPA can be systematically improved being the first 
approximation vi^ithin the so-called adiabatic connection 
fluctuation dissipation (ACFD) framework. 

The RPA is an implicit functional of the density and 
can therefore include non-local correlation effects like e.g. 
the van der Waals interactions. That these are, indeed, 
accurately captured by the RPA has been demonstrated 
in many recent works. ^^"^^ For systems described by 
strong Hubbard-like correlations, the RPA is, however, 
still not fully investigated. A popular test case in this re- 
gard is the dissociation of diatomic molecules with cova- 
Icnt bonds. ^'^ All density functional approximations con- 
structed so far fail in this context if proper spin-symmetry 
is enforced. The total energy in the dissociation limit is 
too high and spurious fractional charges are found at the 
fragments. 

The large error in the total energy has been character- 
ized as static correlation error or so-called fractional spin 
error, studied in detail in the pioneering works of Cohen, 
Mori-Sanchez and Yang.^^"^'^ It has been demonstrated 
that the RPA strongly improves the dissociation limit for 
homoatomic systems such as the H2 molecule. ^^"^^ 

Spurious fractional charges, on the other hand, appear 
only in the dissociation limit of heteroatomic molecules 
and are related to an incorrect behavior of the total en- 
ergy as a function of particle number. The exact en- 
ergy functional exhibits a kink or derivative disconti- 
nuity at integer particle numbers along with a straight 
line behavior between the integers.^® The smooth and 
non-linear behavior of approximate functionals leads to 
charged fragments in the dissociation limit, meaning that 
one electron is too delocalized, i.e., spread over both frag- 



ments. This error is known in the literature as delocaliza- 
tion error or fractional charge error. ^^'^'^'^^ The delocal- 
ization error of the RPA has been studied only for the dis- 
sociation of open-shell and HeJ.'^° It was found to be 
rather severe leading to too low total energies. Whether 
the RPA suffers from fractional charge error in the cases 
of hetcronuclear molecules is presently unknown. 

The exact functional ensures neutral dissociation frag- 
ments by virtue of the KS potential. ^^"'^'^ In the dissocia- 
tion limit the highest occupied orbital of each of the frag- 
ments must be aligned (or degenerate). Consequently, 
the highest occupied molecular orbital (HOMO) of the 
whole system (including both fragments) is a linear com- 
bination of the orbitals of each fragment. To obtain equal 
weights, i.e. integer charge, at the fragments the KS po- 
tential exhibits a sharp step at the bond midpoint, shift- 
ing the energy levels of only one of the fragments. This 
feature is a direct consequence of the derivative disconti- 
nuity in the correlation part of the energy. '^^■'^^ Another 
feature of the exact KS correlation potential of dissoci- 
ated molecules is a peak at the bond midpoint. The 
peak emerges with increasing inter atomic distance, and 
acts to further localize the electrons. 

The RPA correlation potential for atoms was recently 
obtained'^^'^^ and showed a close resemblance to the ex- 
act correlation potential. However, all RPA calculations 
on molecules have so far been carried out using potentials 
originating from other functionals and hence precluding 
a full assessment of the RPA. The aim of this work is 
to provide a more complete analysis of the RPA. To this 
purpose, we have calculated self-consistent RPA poten- 
tials for molecules and investigated the RPA correlation 
potentials for H2 and LiH at different inter-atomic dis- 
tances. These systems allow us to study both the static 
correlation error and the delocalization error. Moreover, 
the LiH is an example where a self-consistent calculation 
is essential. We have also analyzed the RPA energy func- 
tional in terms of fractional charge by studying the RPA 
functional on an extended domain of spin-compensated 
densities allowing for non-integer number of particles. 

We conclude that the RPA potentials exhibit the peak 
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at the bond midpoint but lack the step feature, where 
the latter is related to a missing intra-shell derivative 
discontinuity. The total energy of LiH is, however, still 
largely improved in the dissociation limit as compared 
to, e.g., the exact-exchange (EXX) functional, suggesting 
only a smaller delocalization error. 



II. RPA CORRELATION ENERGY AND 
POTENTIAL 

Within the ACFD framework the exact correlation en- 
ergy is written as'^^"''^ 



Here, we have used the notation (ri,ii) = 1 etc. and 
introduced A(3,2;l) = -iG^(3, 1)G',,(1, 2). The correla- 
tion part of the self-energy Ec in the RPA is given by 



Er, = i- 
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In the appendix we show the expression for ric, i.e., the 
right hand side of Eq. (5), in terms of KS orbitals and 
KS eigenvalues. 
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Tr{«[xAH-XsH]} (1) 



where Xs is the non-interacting KS density response func- 
tion and x\ is the scaled interacting density response 
function. The scaling refers to a system with a linearly 
scaled Coulomb interaction Au(r,r') plus a fictitious po- 
tential which keeps the density fixed as A is changed. The 
parameter A runs between (non-interacting KS case) 
and 1 (fully interacting case). We have used the short 
hand notation Ti fg = J drdr' f{r, r')g{r' , r) for any two- 
point functions / and g. Within TDDFT the function 
reads^^ 



XA = Xs + Xs [Av + /4] XA- 



(2) 



The scaled XC kernel /^c is defined as the functional 
derivative of the scaled XC potential v^^ with respect to 
the density n, evaluated at the ground state density. 

In the RPA f^^ = and thus corresponds to the 
simplest approximation within the ACFD formalism. 
Within the RPA the A-intcgral in Eq. (1) can be evalu- 
ated analytically with the result 
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Diagrammatically Eq. (3) is equal to an infinite summa- 
tion of ring-diagrams. 

The RPA correlation potential Vc can be obtained as 
the functional derivative of Eq. (3) with respect to the 
density. If we let Vs signify the total KS potential, Gg the 
non-interacting KS Green function and Xs = —iGsGs, 
the functional derivative can be obtained via the chain 
rule 
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The result is the well-known linearized Sham-Schliitcr 
(LSS) equation43'44 

/ Xs(l,2)^;c(2)d2= / A(3,2;l)Ee(2,3)d2d3. (5) 



III. FRACTIONAL CHARGE AND SPIN 

The RPA functional produces accurate dissociation en- 
ergies for H2, in contrast to all common density function- 
al which yield a far too high energy due to a spurious 
self-interaction in the H fragments. It is well-known that 
EXX is self-interaction free in the case of a spin-polarized 
H atom. In the dissociation limit of II2, the H atoms are, 
however, not spin-polarized, but rather described by a 
mixture of a spin-up and a spin-down H atom, in which 
case EXX does suffer from self-interaction. To obtain the 
correct dissociation limit a significant correlation contri- 
bution is thus needed. In the following we will show that 
the RPA correlation functional exactly cancels the spuri- 
ous self-interaction in the EXX functional in the dissoci- 
ation limit. The total energy will still not be exact as the 
correlation energy contains a self-correlation term, which 
does not vanish in the dissociation limit. 



A. RPA in the dissociation limit 

The RPA energy is an explicit functional of the KS 
Green function Gg and for spin-compensated systems Gs 
is proportional to the identity matrix in spin-space. In 
the frequency domain a spin-component of Gs reads 

occ uocc 

G,(r, r', c.) = ^ G- (r, r', c.) + G+ (r, r', c.) (8) 

k k 



where we have defined 



G±(r,r',a.) 



V?fc(r)yfc(r^) 
uj- ek±iri 



(9) 



with ipk being a KS spin-orbital and Sk the correspond- 
ing eigenvalue. Consider now a stretched homonuclcar 
diatomic molecule composed of atom A and atom B. For 
large but finite inter-atomic separation R the molecular 
KS orbitals can approximately be written as symmetric 
and antisymmetric linear combination of the atomic KS 
orbitals 
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where (r) is a KS orbital of atom A. This expression 
becomes exact as i? — > 00. It is easy to show that in the 
dissociation hmit 
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where 



Gf — Gf. 

kjto 



uocc 



k=iO 



[Go + G+] . (12) 



Here, Gf contains only states of the isolated atom A. 
The orbital fc = is the special orbital of the highest 
occupied state which has to be considered partially occu- 
pied and partially unoccupied. For a homoatomic system 
the fraction is always 1/2. 

In the case of a covalent bonded heteronuclear diatomic 
molecule a similar analysis can be carried out. The only 
difference is that now only the highest occupied and low- 
est unoccupied arise from a degeneracy of the isolated 
atoms. Due to the lack of symmetry we also have to al- 
low for a more general linear combination of KS orbitals 

1 



(^LUMO(r) ^ -^[V^^o^(r)-v/(2^^^(r)] (13) 

^™(r) = i=[v/(2^^o^(r) + VP^o'^(r)],(14) 

where p G [0,2]. The KS Green functions to be inserted 
in Eq. (11) become 

occ uocc p. 

= EG.T + EG^ + fGo +^G+ (15) 



occ 



uocc 



Go+^G+. (16) 



fe#0 



In the exact system the energy assumes its minimum at 
p = I. This is, however, not guaranteed using an ap- 
proximate functional. In fact, most functionals yield a 
so-called fractional charge error on the atomic fragments 
(i.e. the minimum is found at p ^ 1). 



B. Fractional charge 

In this section we derive an expression for the RPA 
correlation energy in terms of fractionally charged spin- 
compensated systems. To this end, we propose ensembles 
of the following form 
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P 



(i-p)|o)(oi + |(it)(ti + u)(;i) 

for p e [0, 1) 

2-p 



(17) 



(It)(ti + u)(;i) + (p-i)in)(ni(i8) 



for p e (1, 2] 



where |0) refers to the ground state of Nq even number of 
electrons, |ti) referes to a spin-compensated A'o -t- 2-state 
and 11) (It)) a spin-down(up) polarized state with A^o + 1 
electrons. The definition of the non-interacting spin-up 
GJ and spin-down G^ Green functions are 

Gl(rt,r't') = {a\T{4,^rt)^^r't')^}\a) (19) 
Gi{rt,r't') = {a\T{4,HrmHr't'y}\a) (20) 

where a can be any of the states |0), | |), | t) or | H) 
in the KS system. Here, T is the time-ordering operator 
and tpf^{rt) the field operator with the property of adding 
-0''(rt)'l' or removing tjj^{rt) an electron with spin /3 in 
r at time t. Evaluating the Green functions as ensem- 
ble expectation values using Eqs. (17-18) we find Green 
functions exactly of the form as in Eq. (15-16). Thus, 
the ensembles proposed correspond to how an atom in 
the RPA is described in the dissociation limit. 

The RPA correlation energy is usually evaluated using 
the time-ordered KS response function Xs- Using the 
ensemble Green functions we find 



xf(r,r',w) 



~t7TS{u:)p{2~p)\Mr)\^\MrT 
+X?(r,r',w) (21) 



where 

X?(r,r',c.) 



/,(r)/,(r') /,(r)/,(r') 



(22) 



and we have performed a spin summation. The or- 
bitals ipk are now referring to the orbitals of an iso- 
lated atom. The index q = (k, k') is a composite in- 
dex of occupied (fc) and unoccupied (fc') states. The 
special transition q = (0,0) is excluded since it has 
been taken out explicitly (first term in Eq. (21)). The 
functions fq{r) = ipk{v)(pk' (j^) are called excitation func- 
tions and we note that when ipo is occupied it should 
be multiplied by \fpj2 and when it is unoccupied by 
\/{2— p)/2. When calculating the correlation energy the 
interacting RPA response function has to first be evalu- 
ated X\ ~ xf + ^X^'^Xx- The final expression of the 
RPA correlation energy reads 

Er. = -E^l^ f drdv'\^o{r)\'v{r,r')\ipo{v')\' 



pi pOO 1 

+ ^ X y ^^Tr{«[x^(c.)-x?(^)]} (23) 

where x\{uj) is the scaled interacting response function 
calculated using Eq. (22). This equation will allow us 
to investigate the fractional spin error in the following 
section and to obtain numerical results for the fractional 
charge error in Section V. 



C. Fractional spin error 

The fractional spin error is defined as the energy differ- 
ence of a system with A'o -f- 1 electrons with one electron 
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spin polarized (sp) on one hand and the same system but 
fuhy spin compensated (sc) on the other hand. Let us 
spht the RPA interaction energy into the sum of Hartree 
and exchange energy (Hx) and the correlation energy. 
The sum of Hx energy in the sp case reads 



Hx 



No/2 No/2 



(24) 



anti- 



where we introduced the 

symmetrized integral {ij\\kl) = 

J drdr'ip*{r)ip*{r')v(r,r') {2ipk{r)ipi{r^) - ipi{r)ipkir')). 
The second term accounts for all interactions of the 
singly occupied atomic orbital tpo- The expression for 
the spin compensated case reads 



Hx 



pHx 



rfrdr'|(^o(r)|^«(r,r')|^o(r')l'. (25) 



The additional term is a nonzero self-interaction term, 
which is equal to the fractional spin error of EXX. It is 
this term alone that is responsible for the wrong disso- 
ciation limit of EXX for H2. It also introduces a large 
additional error in the dissociation limit of any system, 
on top of the error inherent in the EXX functional when 
describing the fragments. 

We now turn to the RPA correlation energy. The cor- 
relation energy for the sc system is taken from Eq. (23) 
with p = 1. 



El 



1 



drdr'|^o(r)pz;(r,r')|v'o(r')l' 



dX 



duj 
2^ 



TT{v[x\{uj)^xli^)]} (26) 



We now see that the first term exactly cancels the frac- 
tional spin error of EXX. The second term in Eq. (23) is 
identical to the correlation energy obtained in the spin- 
polarized case. Consequently, the RPA functional docs 
not suffer from fractional spin error. There will, how- 
ever, still be a self-correlation due to the second term 
but this is, in general, expected to be smaller. These 
conclusions are confirmed by the numerical results ob- 
tained previously^"'"^^ and here in Sec. V. 



As a reference potential we use the sum of the nuclear 
potential and the Fermi- Amaldi potential, 11-^^ (r), eval- 
uated with the Hartree-Fock density 



.FA 



(r) 



N - 1 
N 



n(r') 



-dr'. 



(28) 



For closed shell systems the derivative of the total en- 
ergy functional with respect to the expansion coefficients, 
bt, is readily evaluated via Eq. (4). For systems with 
fractional charge we evaluated the gradient with the fi- 
nite difference method. The gradient is then fed to the 
Broyden-Fletcher-Goldfarb-Shanno optimization routine 
to find the minimum total energy. The calculations were 
considered converged when the vector norm of the gradi- 
ent was less than 10""^. The algorithm was implemented 
in a local version of the PyQuante module. 

It is a well-known fact that finite basis set OEP cal- 
culations can become numerically unstable if a too large 
potential basis set is used.''''"^^ We carefully chose the po- 
tential basis sets to be balanced to the respective orbital 
basis sets. The proper choice is reflected in the smooth 
potentials that we obtain. For our calculations we used 
the standard orbital basis cc-pVTZ from the Dunning 
family. 

For the potential bases we used even tempered gaus- 
sians. Each set is characterized by three numbers: the 
smallest exponent, amin, the largest exponent, amax, and 
the number of basis functions, N. With this set of pa- 
rameters we construct the exponent as 



(29) 



where i runs from 1 to TV. For the atomic calculations we 
use only s-type functions. For molecular calculations we 
add a set of p-type functions using the same exponents 
and omitting the largest. 

For the one dimensional (ID) calculations we use a 
basis set of cubic splines, which permits us to solve 
Eq. (5) by inverting the static KS response function. 
The spline-basis is described in detail in several previous 
works where it has been used successfully for OEP-type 
of equations. ^'''^"■^^ 



RESULTS 



IV. COMPUTATIONAL DETAILS 

Our goal is to find the local potential that mini- 
mizes the RPA energy functional. For the three dimen- 
sional (3D) calculations we utilize the direct minimiza- 
tion scheme for the optimized effective potential (OEP) 
proposed by Yang and Wu."'^ The potential is expanded 
in a basis plus a reference potential 

v{r)=vo{r)+Y,bt9t{r). (27) 
f 



As a first step we verify our implementation. To this 
end, we compare the total energy and the correlation 
potential with accurate reference data for He, Be, and 
Ne.'^'^ The total energies and ionization potentials (IP) 
are listed in table I. The potential basis sets arc given by 
Omin = 0.01, ttniax = 10, and = 5 for He, amin = 0.01, 
Omax = 10, and A^ = 9 for Ne, and amin = 0.001, amax = 
100, and A^ = 11 for Ne (compare equation (29)). 

The total energies of our new implementation are 
somewhat higher than the reference energies. The dif- 
ferences range from 9 mH to 158 mH. This result is ex- 
pected since the reference calculations were performed 
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this work 


benchmark 


Atom 


IP energy 


IP energy 


He 


-0.907 -2.936 


-0.902 -2.945 


Be 


-0.349 -14.681 


-0.354 -14.754 


Ne 


-0.773 -128.945 


-0.796 -129.103 



RPA correlation potential for Be 



TABLE I: Energies in Hartree. 



RPA correlation potential for He 




10^ 10" 
Distance to nucleus in Bohr 



FIG. 1: The correlation potential of He for RPA^** (dot- 
ted), for RPA in this work (dashed) and the exact correlation 
potential^^ (solid). 



with a virtually complete orbital basis. The differences 
in IP are minor (from 5 mH to 23 mH), thus indicating 
an accurate potential. This conclusion is supported by 
figures 1-3, which show our RPA correlation potentials 
(dashed) of He, Be, and Ne comparing to accurate RPA 
correlation potentials (dotted).'^* The RPA correlation 
potentials shown in the figures are calculated as 



,RPA 



(r) 



RPAr 

-'ks 



riRPA](r)-i;^t Kxx](r). (30) 



In all figures we qualitatively reproduce the reference po- 
tentials. The largest deviation is found close to the nu- 
cleus (the X-axis is on a logarithmic scale). For com- 
parison we also show the exact KS correlation potential 
(solid) and we see that the RPA correlation potential 
closely resembles the exact correlation potential. 

To analyze the RPA molecular correlation poten- 
tial we first investigate a ID model system with 
soft coulomb interactions. The nuclear potentials are 
—Zj^ (x ± i?/2)2 -I- 1, where Z is the nuclear charge and 
i? is the internuclear distance. The electron-electron in- 
teraction is set to 1/ -y/ (x\ — x^)^ + 1. The first system 
consists of two H atoms (Z = 1.2). Figure 4 shows the 
RPA correlation potentials along the bond axis at bond 
distances of 2, 4, and 6 Bohr, with the bond midpoint at 
zero. We notice that the minimum of the correlation po- 
tential is shifted away from the atom, but that the shift 
decreases with increasing interatomic distance. At the 
bond midpoint a positive peak emerges with increasing 
bond distance. Both features are also observed for the 
exact KS correlation potential. ^^■'^^ 




10^ 10" 
Distance to nucleus in Bohr 



FIG. 2: The correlation potential of Be for RPA^* (dotted), 
for RPA in this work (dashed), and the exact correlation 
potential^^ (solid). 



0.15 



RPA correlation potential for Ne 
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FIG. 3: The correlation potential of Ne for RPA^* (dotted), 
for RPA in this work (dashed), and the exact correlation 
potential^^ (solid). 



With our new implementation we are able to investi- 
gate the H2 molecule in all three dimensions and with 
the full coulomb interaction (potential basis: amm = 0.1, 
ttmax = 1.0, and iV = 5). In Figure 5 we show the RPA 
correlation potential (see equation (30)) along the bond 
axis, where again the bond midpoint is at zero. We show 
the potentials for equilibrium distance (1.4 Bohr; solid 
curve), for 5.0 (dotted), and for 10.0 Bohr (dashed). The 
potentials qualitatively resemble the potentials obtained 
from the ID calculations (c.f. 4). Only the peak at the 
bond midpoint is absent for small (1.4) and large (10.0) 
interatomic separation. For large separation the absence 
is easily explained. The orbital and potential basis func- 
tions are simply not diffuse enough to extend to the bond 
midpoint. This results in a vanishing correlation poten- 
tial at the bond midpoint. The situation is different for 
a bond distance of 1.4 Bohr. With a small atomic sepa- 
ration the orbital and potential basis extend to the bond 
midpoint as is evident from the non- vanishing correlation 
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RPA correlation potential of 



ID vl"'"'^ 2 Bohr 
ID 6 Bohr 




2 4 6 8 

Distance to bond midpoint in Bohr 
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0.6 
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RPA correlation potential of LiH 



ID RPA Rnhr 

" v^""'-^ at 6 Bohr 

- - - v^"'"'^ at 8 Bohr 


1 ' 
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-5 5 

Distance to bond midpoint in Bohr 



FIG. 4: The correlation potential along the bond axis with 
the bond midpoint at zero. The system is ID H2 with a soft 
coulomb potential. 



RPA correlation potential of Hj 
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Distance to bond midpoint in Bohr 

FIG. 5: The correlation potential of H2 along the bond axis. 
The bond midpoint is at zero. We show the RPA correlation 
potential for interatomic distance 1.4 (solid), 5.0 (dotted), 
and 10.0 (dashed). 



potential. However, the potential basis functions are not 
compact enough to produce a peak. We have verified this 
by using more compact potential basis functions. With 
this set of basis functions, however, we find unphysical 
oscillations for larger atomic separations. We also placed 
orbital and potential functions at various points along the 
bond axis. We observed a peak in the potential at each 
of the points, which leads us to conclude that these peaks 
are artifacts of the basis rather than genuine features of 
the functional. 

At this point we restrain from showing the dissociation 
curve and rather point the reader to a future publica- 
tion with a detailed discussion of the dissociation curves. 
We would only like to mention that the well-known 
" bump" ^''"^'^ is still present, but somewhat weaker. 

We now turn to the LiH molecule. In figure 6 we show 
the RPA correlation potential (defined in equation (30)) 
of ID LiH (Zli = 3.6, Zh = 1.2). The solid, dotted and 
dashed curves represent the correlation potentials for 2, 
3, and 8 Bohr interatomic distances, respectively. The 
build up of the peak at the bond midpoint is apparent. 



FIG. 6: The correlation potential along the bond with the 
bond midpoint at zero. The system is ID LiH with a soft 
coulomb potential. The Li atom is located at -1, -1.5, and -4 
Bohr, respectively. The H atom is located at 1, 1.5, and 4 
Bohr, respectively. 



However, a step, as is present in the exact KS correlation 
potential, is not observed. 

Figure 7 shows the correlation potential (equa- 
tion (30)) of the three dimensional LiH for bond distances 
3 (solid), 6 (dotted), and 8 Bohr (dashed). The parame- 
ters for the potential basis of the H atom arc amin — 0.01, 
ttmax = 1-0, and N = 2. For the Li atom they read 
ttmin = 0.001, Qfniax = 10, and TV = 3. The same features 
as in the ID case are also found in the results for the 3D 
correlation potentials. In the region of the Li atom (-10 
to 0) the potential qualitatively resembles that of the ID 
system in Figure 6. We see a well with a peak close to 
the nucleus. Also in the region of the H atom (0 to 10) 
figures 6 and 7 show the same structure. A well emerges 
with increasing bond distance. The difference between 
ID and 3D is found only at the bond midpoint. Like in 
the case of H2 a peak emerges only for the ID system 
(figure 6). In contrast, the 3D system (figure 7) exhibits 
a peak at zero only for small and intermediate bond dis- 
tances (3 and 5 Bohr). The explanation for the absence 
of the peak at the bond midpoint for 8 Bohr is the same 
as in the case of H2. 

We further analyze the RPA functional in the context 
of fractional charge. The total energy of H and Li (po- 
tential bases like in LiH) is plotted (figures 8 and 9) as a 
function of the number of electrons, where we allow frac- 
tional values, according to Eq. 23. Both atoms show 
a smooth behavior, thus showing that the RPA (red) 
misses the kink at integer electron numbers. Comparing 
to the exact curves we see, however, a large improvement 
compared to the EXX (blue) functional regarding the to- 
tal energy. The agreement is particularly well at integer 
number of electrons. This explains the good dissociation 
energies for homoatomic systems found in the RPA. We 
can also combine these two figures to analyze the RPA 
energy of LiH in the dissociation limit. In order to do 
this we add the energies for the H atom to the energies 
of the Li atom so that the total number of electrons sums 
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FIG. 7: The correlation potential of LiH along the bond axis. 
The bond midpoint is at zero. The Li atom is located at -L5, 
-2.5, and -4 Bohr, respectively. The H atom is located at 1.5, 
2.5, and 4 Bohr, respectively. 
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FIG. 9: The total energy of Li as a function of number of 
electrons for EXX (blue), RPA (red), and exact (black). 
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FIG. 8: The total energy of H as a function of number of 
electrons for EXX (blue), RPA (red), and exact (black). 
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Number of electrons at the H atom 



FIG. 10: The total energy of dissociated LiH as a function 
of number of electrons at the H atom for EXX (blue) , RPA 
(red), and exact (black). 



up to four. Figure 10 shows the total energy as a function 
of the number of electrons at the H atom. The number 
of electrons at the Li atom will then be four minus the 
X value. The exact functional (black) has a minimum 
at 1.0, because it dissociates LiH into a neutral H atom 
and a neutral Li atom. In Figure 10 we see that EXX 
(blue) and RPA (red) do not dissociate LiH into the neu- 
tral atoms. In both cases there is a surplus of electronic 
charge at the H atom. However, in the case of RPA the 
surplus (0.16 electrons) is much smaller than in the case 
of EXX (0.4 electrons). The large improvement may be 
related to the peak that is present in the RPA correlation 
potential. 

Finally, in table II we collect the total energies in the 
dissociation limit of H2, Li2, and LiH for the exact func- 
tional, EXX, and RPA. Please note that the LiH energies 
of EXX and RPA are taken as the minimum of the re- 
spective curves in figure 10 and not the values at 1.0. 
This is to account for the fractional charge error present 
in both functionals. 



VI. CONCLUSIONS 

We have obtained self-consistent RPA correlation po- 
tentials for diatomic molecules and studied their behavior 
as a function of interatomic distances. At large distances 
the RPA potential correctly exhibits a peak at the bond 
midpoint but misses the step feature. 

We have also analyzed the RPA functional for frac- 
tional charges by evaluating the ensemble averaged KS 
Green function using spin-compensated ensembles with 
non-integer number of particles. This procedure can eas- 
ily be carried over to any functional constructed from the 





exact 


EXX 


RPA 


H2 


-1.000 


-0.713 


-1.035 


Li2 


-14.948 


-14.749 


-14.961 


LiH 


-7.974 


-7.759 


-8.007 



TABLE II: Total energies in Hartree in the dissociation limit. 
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KS Green function. 

The numerical results show that the kink at integer 
number of electrons is missed in the RPA. As a conse- 
quence, we find spurious fractional charges on the dissoci- 
ated fragments. The charges are, however, much smaller 
compared to other functionals. On the other hand, the 
RPA will most likely not be able to describe the so-called 
field counteracting effect in the correlation potential of 
hydrogen chains which has been discussed to have its 
origin in the derivative discontinuity.^'^ 

In summary, we have found that the RPA accomplishes 
the following: 

• The dissociation limit of closed-shell molecules is 
well reproduced in the RPA, due to the explicit 
elimination of the self-interaction term present in 
the EXX functional. 

• RPA separates the charges in bond dissociation by 
virtue of a peak at the bond midpoint, an exact 
feature of the true correlation potential. 

• RPA exhibits only a small fractional charge er- 
ror in the cases of closed-shell covalently bonded 
molecules. 

These findings consolidate the high expectations on the 
RPA that are currently prevalent. There is, however, 
still room for improvement as the discontinuity at odd 
integer particle number is missing. We believe that im- 
provements on the RPA within the ACFD framework will 
help to overcome this shortcoming. 



Appendix A: Derivative of the RPA energy 

In this Appendix we return to Eq. (5) and evaluate 
ric(r) in terms of KS orbitals, tpfc, and KS eigenvalues, 
Efc. The self-energy (Eq. (6)) involves an integration 
over the frequency of the following form 

Sc(^) = I J ^Gsiuj' + uj)vx''''^iu:')v, (Al) 

where we have suppressed the space-coordinates. To per- 
form the integration we write the KS Green function in 
its Lehmann representation 



Gs(r,r',a;) = ^ i^fc(r)^fc(r') 



nk 



1 - 



,(A2) 



where rik is the occupation number. The response func- 
tion in the RPA is given by Eq. (7). This equation can 
easily be rewritten as an eigenvalue problem in terms of 
the matrix 



V, 



qq' 



^qSqq' + {fqWq'), 



(A3) 



where /g(r) = 2y/iU^fq{r), with /g(r) being a KS excita- 
tion function and ujq a KS excitation energy. The square 



root of the eigenvalues corresponds to the true excitation 
energies Zq and the excitation functions are transformed 
according to 



q' 



(A4) 



where U is the matrix diagonalizing V. The interacting 
time-ordered response function in the RPA can then be 
written as 



X«P^(r,r',^) = ^^F,(r)^^,(r') 



1 



U! — Zq + il] UJ + Zq — irj 

and the integral in Eq. (Al) becomes 

Sc(r,r',cj) = ^(^fc(r) J driw(r, ri)F,(ri) 

kq 



(A5) 



2Z„ 



xipk{r') J driv{r',ri)Fq{ri) 

l-Uk , nk 



UJ - Ek - Zq + irj UJ - Ek + Zq - irj 



Next, we evaluate 



duj , 



-i I — I]c(w)Gs(i^)Gs(w). 



(A6) 



(A7) 



After performing standard contour integrations we get in 
total six terms, which after symmetry considerations can 
be reduced to four. In summary, we find 

"c(r) = ^^^(...)</?:(r)¥5p(r) 

ksp q ' 

X J dridr2Fq{ri)v{ri,r2)(pk{r2)f*p{r2) 

X J dridr2ips{ri)tpU^i)v{ri,r2)F*{r2). (A8) 

where the dots signifies the insertion of the following four 
terms 



(1 - nk)npns 



{Ek +Zq- Ep){Ek + Zq- Es) 

^ 2(1 - nk)np{l ~ ris) 

[Ep - Zq- Ek){Ep - Es) 

^ 2nknp{l - Us) 

{Es + Zq- Ek){Ep - Es) 

Ukjl - ng)(l - Up) 

{Es + Zq- Ek){Ep + Zq- Ek) ' 



(A9) 
(AlO) 
(All) 
(A12) 



The expression in Eq. (A8) has, in this work, been im- 
plemented both in the 3D and in the ID case. 
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